1. Field of the Invention
The present invention is directed to a method for motion correction given series of two-dimensional images of a rigid body, wherein the rigid body does not lie entirely in the coverage area of the images and is cut off at an image edge, or at two image edges that do not lie opposite one another.
2. Description of the Prior Art
A number of methods of motion correction for series of time-successive images (scenes) are known from the literature. Such scenes can, for example, be acquired by imaging systems of medical diagnostics such as computed tomography or magnetic resonance tomography, and include movements of an examination subject, for example a patient, during the registration of the scene, which can cause disturbances in the image that are intended to be eliminated by these motion correction methods.
Iterative methods are frequently employed for this purpose, in which all parameters of the motion are simultaneously involved, but these require a high computational outlay.
A Fourier transformation of the image data with a following phase correlation is often applied as an auxiliary technique for the motion determination and correction. For example, European Application 0 368 746 and British Specification 2 187 059 disclose such techniques. German PS 195 00 338 discloses a method for determining the transformation parameters for automatic image registration, wherein the translation vector and rotational angle are determined by one-dimensional correlation of radial projections of selected search image areas from a current image and from a reference image.
The publication of Ohm, J.-R., Digita""e Bildcodierung, Springer-Verlag, 1995, pages 32-34, 42-50, is concerned with the problem of xe2x80x9cedge continuation of image signalsxe2x80x9d, which arises when a filter mask is used that does not coincide with the image coverage area, and thus values exist that lie outside the image coverage region. Various methods are proposed for definition of the signal values at these positions, namely periodic continuation, symmetrical continuation and value-constant continuation.
A further approach is the k-space method for motion correction. This shall be explained in greater detail for the two-dimensional case of a rigid body that lies in the xy-plane. Each movement is then unambiguously identified by the parameters representing translation in the x-direction, translation in the y-direction and rotation around the z-axis. By contrast to the iterative methods, the k-space method decouples the individual parameters from one another. This has the advantage that the sought parameters can be identified more simply by splitting the identification into discrete problems.
Since the Fourier transformation with its properties is the basis of this motion correction method, the relevant definitions and properties of Fourier transformation are first discussed herein, as follows.
Fourier transformation is a linear, mapping rom one mathematical domain to another that is defined as follows in the continuous, two-dimensional case:                                           f            ∼                    ⁡                      (                          k              →                        )                          =                              ∫                                          X                →                                            ϵ                                  R                  2                                                              ⁢                                                    f                ⁡                                  (                                      x                    →                                    )                                            ·                              ⅇ                                  2                  ⁢                  m                  ⁢                                      xe2x80x83                                    ⁢                                      k                    →                                    ⁢                                      x                    →                                                                        ⁢                          xe2x80x83                        ⁢                          ⅆ                              x                →                                                                        (        1        )                                                      f            ~                    ⁡                      (                          X              →                        )                          =                              ∫                                          X                →                                            ϵ                                  R                  2                                                              ⁢                                                    f                ⁡                                  (                                      x                    →                                    )                                            ·                              ⅇ                                                      -                    2                                    ⁢                  m                  ⁢                                      xe2x80x83                                    ⁢                                      k                    →                                    ⁢                                      x                    →                                                                        ⁢                          xe2x80x83                        ⁢                          ⅆ                              x                →                                                                        (        2        )            
The target domain of the Fourier transformation is called k-space and the source domain is called the spatial domain. When f(x) is real, then the k-space becomes point-symmetrical relative to the zero point upon application of complex conjugation:                                           ∫                                          x                →                                            ϵ                ⁢                                  xe2x80x83                                ⁢                                  R                  2                                                              ⁢                      f            ⁢                          xe2x80x83                        ⁢                                          (                                  x                  →                                )                            ·                              (                                  ⅇ                                      2                    ⁢                    m                    ⁢                                          xe2x80x83                                        ⁢                                          (                                                                        -                                                      k                            →                                                                          ⁢                                                  xe2x80x83                                                ⁢                                                  x                          →                                                                                                                    )                                      *                          xe2x80x83                        ⁢                          ⅆ                              x                →                                                    =                                                            ∫                                                      x                    _                                                        ϵ                    ⁢                                          xe2x80x83                                        ⁢                                          R                      2                                                                                  ⁢                              f                ⁢                                  xe2x80x83                                ⁢                                                      (                                          x                      →                                        )                                    ·                                      ⅇ                                          2                      ⁢                      m                      ⁢                                              xe2x80x83                                            ⁢                                              k                        →                                            ⁢                                              xe2x80x83                                            ⁢                                              x                        →                                                                                            ⁢                                  xe2x80x83                                ⁢                                  ⅆ                                      x                    →                                                                        ⇔                                          f                ~                            ⁢                              xe2x80x83                            ⁢                              (                                  k                  ~                                )                                              =                                    f              ~                        ⁢                          xe2x80x83                        ⁢                          (                              -                                  k                  →                                            )                        *                                              (        3        )            
When x is shifted by an arbitrary vector, then the function f(x+a) can be calculated using the Fourier transform of f(x):                                                                         f                ⁡                                  (                                                            X                      →                                        +                                          a                      →                                                        )                                            =                                                ∫                                                            k                      →                                        ∈                                          R                      2                                                                      ⁢                                                                                                    f                        →                                            ⁡                                              (                                                  k                          →                                                )                                                              ·                                          ⅇ                                                                        -                          2                                                ⁢                        m                        ⁢                                                  xe2x80x83                                                ⁢                                                  k                          →                                                ⁢                                                  xe2x80x83                                                ⁢                                                  (                                                                                    x                              →                                                        =                                                          a                              →                                                                                )                                                                                                      ⁢                                      xe2x80x83                                    ⁢                                      ⅆ                                          k                      →                                                                                                                                              =                                                ∫                                                            k                      →                                        ∈                                          xe2x80x83                                        ⁢                                          R                      2                                                                      ⁢                                                                                                    f                        →                                            ⁡                                              (                                                  k                          →                                                )                                                              ·                                                                  ⅇ                                                                              -                            2                                                    ⁢                          m                          ⁢                                                      xe2x80x83                                                    ⁢                                                      k                            →                                                    ⁢                                                      a                            →                                                                                                                      ⏟                                                  ϕ                          ⁡                                                      (                                                          k                              →                                                        )                                                                                                                ·                                          ⅇ                                                                        -                          2                                                ⁢                        m                        ⁢                                                  xe2x80x83                                                ⁢                                                  k                          →                                                ⁢                                                  x                          →                                                                                                      ⁢                                      ⅆ                                          k                      →                                                                                                                              (        4        )            
A shift of the image in the spatial domain is thus expressed in an additional phase xcfx86({right arrow over (k)}) of the Fourier transform. A rotation in the location space by an angle xcex8 expressed by multiplication with a rotation matrix (Rxcex8), yields:                                                                                           f                  ⁡                                      (                                                                  R                        ⁡                                                  (                          θ                          )                                                                    ·                                              X                        →                                                              )                                                  =                                                      ∫                                                                  k                        →                                            ∈                                              xe2x80x83                                            ⁢                                              R                        2                                                                              ⁢                                                                                                              f                          →                                                ⁡                                                  (                                                      k                            →                                                    )                                                                    ·                                              ⅇ                                                                              -                            2                                                    ⁢                          m                          ⁢                                                      xe2x80x83                                                    ⁢                                                      k                            →                                                    ⁢                                                      R                            ⁡                                                          (                              θ                              )                                                                                ⁢                                                      x                            →                                                                                                                ⁢                                          xe2x80x83                                        ⁢                                          ⅆ                                              k                        →                                                                                                                                                                    =                                                      ∫                                                                  k                        →                                            ∈                                              xe2x80x83                                            ⁢                                              R                        2                                                                              ⁢                                                                                                              f                          →                                                ⁡                                                  (                                                      k                            →                                                    )                                                                    ·                                              ⅇ                                                                              -                            2                                                    ⁢                          m                          ⁢                                                      xe2x80x83                                                    ⁢                                                                                    R                                                              -                                1                                                                                      ⁡                                                          (                              θ                              )                                                                                ⁢                                                      k                            →                                                    ⁢                                                      x                            →                                                                                              ·                                              ⅇ                                                                              -                            2                                                    ⁢                          m                          ⁢                                                      xe2x80x83                                                    ⁢                                                      k                            →                                                    ⁢                                                      x                            →                                                                                                                ⁢                                          ⅆ                                              k                        →                                                                                                                                ⁢                  
                ⁢                              with            ⁢                          xe2x80x83                        ⁢                          R              ⁡                              (                θ                )                                              =                      (                                                                                cos                    ⁡                                          (                      θ                      )                                                                                                            -                                          sin                      ⁡                                              (                        θ                        )                                                                                                                                                              sin                    ⁡                                          (                      θ                      )                                                                                                            cos                    ⁡                                          (                      θ                      )                                                                                            )                                              (        5        )            
and with the substitution {right arrow over (k)}xe2x80x2=Rxe2x88x921(xcex8){right arrow over (k)}=R(xe2x88x92xcex8){right arrow over (k)}:                                                         ⅆ                              k                z                xe2x80x2                                                    ⅆ                              k                →                                              =                      (                                                                                cos                    ⁡                                          (                                              -                        θ                                            )                                                                                                            sin                    ⁡                                          (                                              -                        θ                                            )                                                                                                                                        -                                          sin                      ⁡                                              (                                                  -                          θ                                                )                                                                                                                                  cos                    ⁡                                          (                                              -                        θ                                            )                                                                                            )                          ⁢                  
                ⁢                              ⅆ                          k              y              xe2x80x2                                            ⅆ                          k              →                                      =                                            (                                                                                          cos                      ⁡                                              (                                                  -                          θ                                                )                                                                                                                        -                                              sin                        ⁡                                                  (                                                      -                            θ                                                    )                                                                                                                                                                                sin                      ⁡                                              (                                                  -                          θ                                                )                                                                                                                        cos                      ⁡                                              (                                                  -                          θ                                                )                                                                                                        )                        ⁢                          
                        ⇒                          ⅆ                                                k                  →                                xe2x80x2                                              :=                      xe2x80x83                    ⁢                                                    ⅆ                                  k                  x                  xe2x80x2                                            ·                              ⅆ                                  k                  y                  xe2x80x2                                                      =                                                            (                                                                                                              cos                          ⁡                                                      (                                                          -                              θ                                                        )                                                                                                                                                sin                          ⁡                                                      (                                                          -                              θ                                                        )                                                                                                                                                                                        -                                                      sin                            ⁡                                                          (                                                              -                                θ                                                            )                                                                                                                                                                            cos                          ⁡                                                      (                                                          -                              θ                                                        )                                                                                                                                )                                ·                                  (                                                                                                              cos                          ⁡                                                      (                                                          -                              θ                                                        )                                                                                                                                                -                                                      sin                            ⁡                                                          (                                                              -                                ϑ                                                            )                                                                                                                                                                                                                    sin                          ⁡                                                      (                                                          -                              θ                                                        )                                                                                                                                                cos                          ⁡                                                      (                                                          -                              θ                                                        )                                                                                                                                )                                ·                                  ⅆ                                      k                    →                                                              ⁢                              
                            ⇒                              f                (                                                                            R                      ⁡                                              (                        θ                        )                                                              ·                                          X                      →                                                        =                                                            ∫                                                                        k                          →                                                ∈                                                  xe2x80x83                                                ⁢                                                  R                          2                                                                                      ⁢                                                                                                                        f                            ∼                                                    ⁡                                                      (                                                                                          R                                ⁡                                                                  (                                  θ                                  )                                                                                            ⁢                                                                                                k                                  →                                                                xe2x80x2                                                                                      )                                                                          ·                                                  ⅇ                                                                                    -                              2                                                        ⁢                            m                            ⁢                                                          xe2x80x83                                                        ⁢                                                          k                              →                                                        ⁢                                                          x                              →                                                                                                                          ⁢                                              ⅆ                                                                              k                            →                                                    xe2x80x2                                                                                                                                                                            (        6        )            
A rotation of the location space by the angle xcex8 effects the same rotation of the Fourier transform.
Overall, thus, for an arbitrary motion:
Given the shift of two functions f and g, Rxe2x86x92R with f(x)xe2x88x92g(x+a),
the following is valid with respect to its Fourier transform:
{tilde over (ƒ)}(k)={tilde over (g)}(k)xc2x7exe2x88x922mka xe2x80x83xe2x80x83(7) 
{tilde over (ƒ)}(k)xc2x7{tilde over (g)}(k)*=|{tilde over (g)}(k)|2xc2x7exe2x88x922mka xe2x80x83xe2x80x83(8) 
                                                                        f                ∼                            ⁡                              (                k                )                                      ·                                                            g                  ∼                                ⁡                                  (                  k                  )                                            *                                                          "LeftBracketingBar"                                                g                  ∼                                ⁡                                  (                  k                  )                                            "RightBracketingBar"                        2                          =                  ⅇ                                    -              2                        ⁢            mka                                              (        9        )            
If the motion parameters of a rigid body in a plane that is arbitrarily shifted and rotated in two successive images is to be found, then this problem can be reduced to a number of one-dimensional observations (shifts per two functions Rxe2x86x92R). The prerequisite for the method, however, is that the rigid body does not extend beyond the coverage area of the images. For example, the movement of a rigid body in the xy-plane of a three dimensional Cartesian coordinate system is unambiguously characterized by the three parameters of translation in the x-direction, translation in the y-direction and rotation around the z-axis.
Determination of the Rotation Angle:
According to equation (6), the rotation of an image in the spatial domain also causes the rotation of the Fourier transform. A translation, however, only acts on the phase of the Fourier transform. In order to decouple the rotation from the translation, the magnitude of the Fourier transform is thus considered. Each rotation here is a rotation around the origin. The sum over an annulus around the origin is invariant to the rotation. When an annulus is summed along its radii, then a one-dimensional function g(xcfx86) is obtained. Two images rotated relative to one another thus supply two oppositely shifted functions g1(xcfx86) and g2(xcfx86). Since the image data in the location space are real, the k-space is symmetrical. Summation over respective halves of the annuluses thus suffices. According to equation (9), the shift can now be identified. When the shift of the two functions g1(xcfx86) and g2(xcfx86) is known, then the rotational angle has been found.
The rotational angle of a rigid body presented in two successive images is thus known and can be correspondingly corrected.
Subsequently, the translation parameters to be considered are determined, these being relative to the rotation and being independent of one another. As long as the prerequisite of the movement of a rigid body is met, i.e. as long as the subject does not leave the image area due to the translation, the translations along the x-axis and the y-axis can be separated from one another. This occurs by projections of the body onto the coordinate axes. These projections are respectively two lines shifted relative to one another. The determination of the translation parameters is thus reduced to the same sub-problem as the determination of the rotational angle, namely the problem of two functions Rxe2x86x92R shifted relative to one another. The solution of this problem also ensues according to equation (9).
The above, known motion correction method is based on a rigid body that moves and/or is rotated from one image to another but does not change in shape. This prerequisite, however, is not satisfied when only a portion of an inherently rigid body lie in the coverage area of the image and other portions extend beyond the coverage area at an image edge. Because the rigid body is moved and/or rotated, parts of the body disappear from the coverage area of the images or new parts are added thereto. The presented body thus changes in shape, and a rigid body can longer be assumed in the observation of the images.
The prerequisites for the initially described method thus are not satisfied. This only occurs, however, in the edge zones of the images, and is not distributed over the entire k-space by the Fourier transformation. As described above, thus, the known method is not suitable under this circumstance and more complicated methods must be used.
An object of the present invention is to provide an efficient motion correction method that can also be used for a series of successive images wherein the body to be imaged does not lie completely in the coverage area of the images and is therefore cut off at the image edge.
The above object is achieved in a first embodiment of the inventive method, wherein a series of two-dimensional images of a rigid body is obtained, the images respectively having image edges defining a coverage area and wherein the rigid body does not lie completely in the coverage area and is cut off by one of the image edges, including the steps of selecting first and second chronologically successive images in the series, between which the motion correction is to implemented, and respectively mirroring the first and second images at the image edge which cuts off the rigid body, thereby producing two mirror images, attaching the first image and the mirror image of the second image and attaching the second image and the mirror image of the first image, at the image edge which cuts off the rigid body, thereby producing two new images that show a rigid body, to a good approximation, lying completely within these new images, determining motion parameters for time-successive images of the rigid body from the two new images according to a motion correction method, and correcting the first and second images using these motion parameters.
The above object is also achieved in a second embodiment of the inventive method, for motion correction in a series of two-dimensional images of a rigid body, the images respectively having image edges defining a coverage area and image edges including two image edges that do not lie opposite each other, and wherein the rigid body does not lie completely in the coverage area of the images and is cut off by the aforementioned two image edges that do not lie opposite each other, including the steps of selecting first and second time-successive images in the series, between which the motion correction is to be implemented, and respectively mirroring each of said two image edges and a point of intersection between these two image edges, thereby producing three mirror images for each of said first and second images, attaching the three mirror images acquired from the second image to the first image respectively at the two image edges, and attaching the three mirror images acquired from the first image to the second image respectively at the two image edges, thereby producing two new images that show a rigid body, to a good approximation, lying completely within these new images, determining motion parameters for time-successive images of a rigid body from the two new images according to a motion correction method, and correcting the first and second images using these motion parameters.
The newly arising images in step b) in each of the above embodiments of the inventive method, between which the further motion correction is to be implemented, exhibit a rigid body that completely lies in the coverage area of these images in a good approximation for small rotations and small translations.
According to the invention, a correction method thus is made available that functions in the manner of methods operating on the assumption of a rigid body that lies completely in the coverage area of the time-successive images. As discussed above, such methods are known and require significantly less calculating outlay compared to calculating methods for which the above premise does not apply, even if the necessary mirrorings and the images that are now two or four times as large require additional calculating outlay.
The series of successive images can be of tomograms as well as of projections of registered volumes.
Given the assumption that the rigid body can be projected into three orthogonal planes, whereby the projections either completely contain the image to be displayed or cut this off at an image edge, or at two image edges not lying opposite one another, the inventive method can also be expanded to three-dimensional bodies. The problem is thus resolved into three mutually independent sub-problems that are then solved by the motion correction methods according to the above embodiments.
The expansion of the motion correction to three dimensions initially increases the number of parameter to be identified from 3 to 6. A body can now be rotated around all three coordinates axes and shifted along all three axes. As in the two-dimensional case, the rotation parameters are first identified and the rotations are correspondingly corrected.
In order to decouple the rotations from the translations, the magnitude of the Fourier transform is again considered. In the two-dimensional case, the sum over an annulus whose center lies on the zero frequency in the k-space is invariant relative to rotations. In the three-dimensional case, this invariant is the sum over each hollow sphere whose mid-point lies on the zero frequency. In this case as well, the projection along the radii is possible, since each rotation after the Fourier transformation and amount formation occurs around the origin. This projection generates two discrete functions on two oppositely rotated spherical surfaces. Given rotation of the spherical surfaces relative to one another, three degrees of freedom are possible that, for example can be described by the three Euler angles known from classical mechanics.
The separation of the parameters from one another is problematical. This cannot be achieved by a simple projection as in the case of the translations. The parameters cannot be separated in this way. Given observation of three annuluses around the coordinate axes, thus, each annulus is also dependent on the rotations around the other axes. Mathematically, this problem is based on the fact that the two three-dimensional rotational matrices are generally not commutative. A transformation into a domain that produces commutativity would be necessary for a parameter separation.
The separation of the parameters, however, nonetheless can be achieved in many applications, however, it is necessary for this purpose to limit the allowable rotational angles. The advantage of the k-space techniques of also being able to recognize large rotations is lost in the three-dimensional case.
For small rotational angles, the rotations around the three axes are considered a good approximation, as being independent of one another. A volume data set is thus projected into the three orthogonal planes: xy-plane, yz-plane and xz-plane. The rotational angles are subsequently determined as described at the outset on the three two-dimensional data sets resulting therefrom. The error that is committed by this approximative solution lies in the 0.1% range given, for example, rotations by 2xc2x0 and thus can be left out of consideration.
When the rotational angles have thus been found and the data sets have been correspondingly corrected, then the determination of the translation parameters follows, as in the two-dimensional case as well. The translation parameters, as in the two-dimensional case, are determined decoupled arid independently from one another by projections of the data sets onto the coordinate axes.
The motion correction for the three-dimensional thus also is been reduced to the shift of two respective functions Rxe2x86x92R and can be solved with the initially described method.
The method of the invention can be especially advantageously utilized in 3-D image registration methods such as computed tomography or magnetic resonance tomography. The reduced calculating outlay is important particularly for the medical application of the method, wherein short image registration and processing times are especially desirable.
In such medical examinations, examinations of the head as well as of the brain represent preferred areas of application of the inventive method because the prerequisite for the method, that the body to be displayed is only cut off at one image edge, or at two image edges not lying opposite one another, can be unproblematically met. The neck region given presentations of the head, or the spinal marrow given presentations of the brain, thus usually leave the coverage area at an image edge.
According to experimental examinations for two time-successive images of a rigid body to be corrected relative one another, the limits of the inventive method lie at approximately 20% of the expanse of the body in the motion direction for the translational motion and at approximately 5xc2x0 for the rotation. Given scans that exceed these limits, the image quality after the correction decreases noticeably. For the majority of medical applications, however, these limits do not represent a practical limitation of the inventive method.